#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <map>
#include <algorithm>
#include <fstream>
#include <time.h>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/tools/Pruner.hh"

#include "Njettiness.hh"

#include "CmdLine.hh"
#include "TFile.h"
#include "TChain.h"
#include "TProfile.h"

#include "TH1D.h"
#include "TH2D.h"

#include "MyAnalysis.h"

void MyAnalysis::hist_init()
{
//initialize histograms
  
  hist["nevt"]=new TH1D("nevt", "nevt", 10, 0, 10);
  hist["nevt_unweighted"]=new TH1D("nevt_unweighted", "nevt_unweighted",
				   10, 0, 10);

  hist["met"]=new TH1D("MEt", "#slash{E}_{T}", 50, 0, 500);

  hist["jetmass"]=new TH1D("jetmass", "jetmass", mass/5, 0, mass*2);
  hist["jetmass1"]=new TH1D("jetmass1", "jetmass1", mass/5, 0, mass*2);
  hist["jetmass2"]=new TH1D("jetmass2", "jetmass2", mass/5, 0, mass*2);
  hist["jetmass3"]=new TH1D("jetmass3", "jetmass3", mass/5, 0, mass*2);

  hist["jetpt"]=new TH1D("jetpt", "jetpt", mass/5, 0, mass*2);
  hist["jetpt1"]=new TH1D("jetpt1", "jetpt1", mass/5, 0, mass*2);
  hist["jetpt2"]=new TH1D("jetpt2", "jetpt2", mass/5, 0, mass*2);
  hist["jetpt3"]=new TH1D("jetpt3", "jetpt3", mass/5, 0, mass*2);


  hist["njet"]=new TH1D("njet", "njet", 40, 0, 40);

  hist["Mt"]=new TH1D("Mt", "Mt", mass*2/5, 0, mass*4);
  hist["Ht"]=new TH1D("Ht", "Ht", mass*2/5, 0, mass*4);
  hist2["Mt_Ht"]=new TH2D("Mt_Ht", "Mt_Ht", mass*2/5, 0, mass*4, mass*2/5, 0, mass*4);
  hist2["Mt_radius"]=new TH2D("Mt_radius", "M_{J} vs. jet radius", 6, 0.6, 1.8, 1000, 0, 2000);

  hist["jet1_tau1"]=new TH1D("jet1_tau1", "1^{st} jet #tau_{1}", 50, 0, 1);
  hist["jet1_tau2"]=new TH1D("jet1_tau2", "1^{st} jet #tau_{2}", 50, 0, 1);
  hist["jet1_tau3"]=new TH1D("jet1_tau3", "1^{st} jet #tau_{3}", 50, 0, 1);

  hist["jet1_tau21"]=new TH1D("jet1_tau21", "1^{st} jet #tau_{2}/#tau_{1}", 50, 0, 1);
  hist["jet1_tau32"]=new TH1D("jet1_tau32", "1^{st} jet #tau_{3}/#tau_{2}", 50, 0, 1);
  hist["jet1_tau31"]=new TH1D("jet1_tau31", "1^{st} jet #tau_{3}/#tau_{1}", 50, 0, 1);

  hist2["jet1_tau3_tau2"]=new TH2D("jet1_tau3_tau2", 
				  "1^{st} jet #tau_{3} vs #tau_{2}", 
				  50, 0, 1, 50, 0, 1);

  hist["jet2_tau1"]=new TH1D("jet2_tau1", "2^{nd} jet #tau_{1}", 50, 0, 1);
  hist["jet2_tau2"]=new TH1D("jet2_tau2", "2^{nd} jet #tau_{2}", 50, 0, 1);
  hist["jet2_tau3"]=new TH1D("jet2_tau3", "2^{nd} jet #tau_{3}", 50, 0, 1);

  hist["jet2_tau21"]=new TH1D("jet2_tau21", "2^{nd} jet #tau_{2}/#tau_{1}", 50, 0, 1);
  hist["jet2_tau32"]=new TH1D("jet2_tau32", "2^{nd} jet #tau_{3}/#tau_{2}", 50, 0, 1);
  hist["jet2_tau31"]=new TH1D("jet2_tau31", "2^{nd} jet #tau_{3}/#tau_{1}", 50, 0, 1);

  hist2["jet2_tau3_tau2"]=new TH2D("jet2_tau3_tau2", 
				   "2^{nd} jet #tau_{3} vs #tau_{2}", 
				   50, 0, 1, 50, 0, 1);
  

  //N-jetiness
  /*
  hist["glu_tau1"]=new TH1D("glu_tau1", "gluino jet #tau_{1}", 50, 0, 1);
  hist["glu_tau2"]=new TH1D("glu_tau2", "gluino jet #tau_{2}", 50, 0, 1);
  hist["glu_tau3"]=new TH1D("glu_tau3", "gluino jet #tau_{3}", 50, 0, 1);
  
  hist["glu_tau21"]=new TH1D("glu_tau21", "gluino jet #tau_{2}/#tau_{1}", 50, 0, 1);
  hist["glu_tau32"]=new TH1D("glu_tau32", "gluino jet #tau_{3}/#tau_{2}", 50, 0, 1);
  hist["glu_tau31"]=new TH1D("glu_tau31", "gluino jet #tau_{3}/#tau_{1}", 50, 0, 1);
  */

  //jettiness in bins of p_T
  
  hist["jet1_tau31_ptbin1"]=new TH1D("jet1_tau31_ptbin1", "First Jet #tau_{3}/#tau_{1} in pt bin1", 50,0,1);
  hist["jet1_tau31_ptbin2"]=new TH1D("jet1_tau31_ptbin2", "First Jet #tau_{3}/#tau_{1} in pt bin2", 50,0,1);
  hist["jet1_tau31_ptbin3"]=new TH1D("jet1_tau31_ptbin3", "First Jet #tau_{3}/#tau_{1} in pt bin3", 50,0,1);
  hist["jet1_tau31_ptbin4"]=new TH1D("jet1_tau31_ptbin4", "First Jet #tau_{3}/#tau_{1} in pt bin4", 50,0,1);
  hist["jet1_tau31_ptbin5"]=new TH1D("jet1_tau31_ptbin5", "First Jet #tau_{3}/#tau_{1} in pt bin5", 50,0,1);


  hist["jet1_tau21_ptbin1"]=new TH1D("jet1_tau21_ptbin1", "First Jet #tau_{2}/#tau_{1} in pt bin1", 50,0,1);
  hist["jet1_tau21_ptbin2"]=new TH1D("jet1_tau21_ptbin2", "First Jet #tau_{2}/#tau_{1} in pt bin2", 50,0,1);
  hist["jet1_tau21_ptbin3"]=new TH1D("jet1_tau21_ptbin3", "First Jet #tau_{2}/#tau_{1} in pt bin3", 50,0,1);
  hist["jet1_tau21_ptbin4"]=new TH1D("jet1_tau21_ptbin4", "First Jet #tau_{2}/#tau_{1} in pt bin4", 50,0,1);
  hist["jet1_tau21_ptbin5"]=new TH1D("jet1_tau21_ptbin5", "First Jet #tau_{2}/#tau_{1} in pt bin5", 50,0,1);


  hist["jet1_tau31_mbin1"]=new TH1D("jet1_tau31_mbin1", "First Jet #tau_{3}/#tau_{1} in mass bin1", 50,0,1);
  hist["jet1_tau31_mbin2"]=new TH1D("jet1_tau31_mbin2", "First Jet #tau_{3}/#tau_{1} in mass bin2", 50,0,1);
  hist["jet1_tau31_mbin3"]=new TH1D("jet1_tau31_mbin3", "First Jet #tau_{3}/#tau_{1} in mass bin3", 50,0,1);
  hist["jet1_tau31_mbin4"]=new TH1D("jet1_tau31_mbin4", "First Jet #tau_{3}/#tau_{1} in mass bin4", 50,0,1);

  hist["jet1_tau21_mbin1"]=new TH1D("jet1_tau21_mbin1", "First Jet #tau_{2}/#tau_{1} in mass bin1", 50,0,1);
  hist["jet1_tau21_mbin2"]=new TH1D("jet1_tau21_mbin2", "First Jet #tau_{2}/#tau_{1} in mass bin2", 50,0,1);
  hist["jet1_tau21_mbin3"]=new TH1D("jet1_tau21_mbin3", "First Jet #tau_{2}/#tau_{1} in mass bin3", 50,0,1);
  hist["jet1_tau21_mbin4"]=new TH1D("jet1_tau21_mbin4", "First Jet #tau_{2}/#tau_{1} in mass bin4", 50,0,1);



  //initialize mass plot as a function of radius
  hist["jetmass_0.6"]=new TH1D("jetmass_0.6", 
			       "Jet mass with radius 0.6",
			       1000, 0, 2000);
  hist["jetmass_0.8"]=new TH1D("jetmass_0.8", 
			       "Jet mass with radius 0.8",
			       1000, 0, 2000);
  hist["jetmass_1.0"]=new TH1D("jetmass_1.0", 
			       "Jet mass with radius 1.0",
			       1000, 0, 2000);
  hist["jetmass_1.2"]=new TH1D("jetmass_1.2", 
			       "Jet mass with radius 1.2",
			       1000, 0, 2000);
  hist["jetmass_1.4"]=new TH1D("jetmass_1.4", 
			       "Jet mass with radius 1.4",
			       1000, 0, 2000);
  hist["jetmass_1.6"]=new TH1D("jetmass_1.6", 
			       "Jet mass with radius 1.6",
			       1000, 0, 2000);
  mass_profile=new TProfile("jetmass_radius",
			    "Jet mass vs. radius",
			    8, 0.4, 2.0);
  
  jetmass_plots.push_back(hist["jetmass_0.6"]);
  jetmass_plots.push_back(hist["jetmass_0.8"]);
  jetmass_plots.push_back(hist["jetmass_1.0"]);
  jetmass_plots.push_back(hist["jetmass_1.2"]);
  jetmass_plots.push_back(hist["jetmass_1.4"]);
  jetmass_plots.push_back(hist["jetmass_1.6"]);
  

}

void MyAnalysis::hist_init2()
{
  hist2["jet1_mass_tau21"]=new TH2F("jet1_mass_tau21",
				    "jet1_mass_tau21",
				    50, 0, 1,
				    50, 0, 500);

  hist2["jet1_mass_tau31"]=new TH2F("jet1_mass_tau31",
				    "jet1_mass_tau31",
				    50, 0, 1,
				    50, 0, 500);

  hist2["jet1_mass_tau32"]=new TH2F("jet1_mass_tau32",
				    "jet1_mass_tau32",
				    50, 0, 1,
				    50, 0, 500);


  hist2["jet1_pt_tau21"]=new TH2F("jet1_pt_tau21",
				  "jet1_pt_tau21",
				  50, 0, 1,
				  100, 0, 2000);

  hist2["jet1_pt_tau31"]=new TH2F("jet1_pt_tau31",
				  "jet1_pt_tau31",
				  50, 0, 1,
				  100, 0, 2000);

  hist2["jet1_pt_tau32"]=new TH2F("jet1_pt_tau32",
				  "jet1_pt_tau32",
				  50, 0, 1,
				  100, 0, 2000);

  hist2["jet2_pt_tau31"]=new TH2F("jet2_pt_tau31",
				  "jet2_pt_tau31",
				  50, 0, 1,
				  100, 0, 2000);

  hist2["jet2_pt_tau21"]=new TH2F("jet2_pt_tau21",
				  "jet2_pt_tau21",
				  50, 0, 1,
				  100, 0, 2000);

  hist2["jet2_pt_tau32"]=new TH2F("jet2_pt_tau32",
				  "jet2_pt_tau32",
				  50, 0, 1,
				  100, 0, 2000);


  hist2["Mt_tau21"]=new TH2F("Mt_tau21",
			     "Mt_tau21",
			     50, 0, 1,
			     100, 0, 1000);
  
  hist2["Mt_tau31"]=new TH2F("Mt_tau31",
			     "Mt_tau31",
			     50, 0, 1,
			     100, 0, 1000);
  
  hist2["Mt_tau32"]=new TH2F("Mt_tau32",
			     "Mt_tau32",
			     50, 0, 1,
			     100, 0, 1000);


  hist2["Mt_m_tau21"]=new TH2F("Mt_m_tau21",
			     "Mt vs tau21 for the most massive jet",
			     50, 0, 1,
			     100, 0, 1000);
  
  hist2["Mt_m_tau31"]=new TH2F("Mt_m_tau31",
			     "Mt vs tau31 for the most massive jet",
			     50, 0, 1,
			     100, 0, 1000);

  hist["Mjjj_100"]=new TH1D("Mjjj_100", "m_{jjj} after 100 diagonal cut", 150, 0, 1500);  
  hist["Mjjj_200"]=new TH1D("Mjjj_200", "m_{jjj} after 200 diagonal cut", 150, 0, 1500);
  hist["Mjjj_300"]=new TH1D("Mjjj_300", "m_{jjj} after 300 diagonal cut", 150, 0, 1500);
  hist["Mjjj_400"]=new TH1D("Mjjj_400", "m_{jjj} after 400 diagonal cut", 150, 0, 1500);
  hist["Mjjj_500"]=new TH1D("Mjjj_500", "m_{jjj} after 500 diagonal cut", 150, 0, 1500);

  hist2["Mjjj_sumpt"]=new TH2D("Mjjj_sumpt", "m_{jjj} vs. #sum p_{T}", 
			       150, 0, 1500, 150, 0, 1500);

  hist2["Mjjj_sumpt_max"]=new TH2D("Mjjj_sumpt_max", 
				   "m_{jjj} vs. #sum p_{T} for highest m_{jjj}-#sum p_{T} triplet",
				   150, 0, 1500, 150, 0, 1500);


  for(int i=0; i<n_display; i++)
    {
      int evt_entry=i;
      key_entry.push_back(evt_entry);

      stringstream sout;
      sout<<"display_"<<i;
      fileout->mkdir(sout.str().c_str())->cd();
      
      sout.str("");
      sout<<"dis_npar_"<<evt_entry;
      key_E.push_back(sout.str());
      hist2[sout.str()]=new TH2D(sout.str().c_str(),
				 sout.str().c_str(),
				 50, -3, 3,
				 50, 0, 3.142*2);
      hist2[sout.str()]->SetOption("lego");
	
      sout.str("");
      sout<<"dis_Pt_"<<evt_entry;
      key_Pt.push_back(sout.str());
      hist2[sout.str()]=new TH2D(sout.str().c_str(),
				 sout.str().c_str(),
				 50, -3, 3,
				 50, 0, 3.142*2);
      hist2[sout.str()]->SetOption("lego"); 

      //grab the first 4 jets
      key_jet.push_back(vector<string>());
      for(int j=0; j<4; j++)
	{
	  sout.str("");
	  sout<<"Pt_"<<i<<"_jet"<<j+1;
	  key_jet[i].push_back(sout.str());
	  hist2[sout.str()]=new TH2D(sout.str().c_str(),
				     sout.str().c_str(),
				     50, -3, 3,
				     50, 0, 3.142*2);
	  hist2[sout.str()]->SetOption("lego");
	  hist2[sout.str()]->SetLineColor(j+2);
	}

      if(b_gluino)
	{
	  sout.str("");
	  sout<<"glu1_"<<i;
	  hist2[sout.str()]=new TH2D(sout.str().c_str(),
				     sout.str().c_str(),
				     50, -3, 3,
				     50, 0, 3.142*2);
	  hist2[sout.str()]->SetOption("lego");
	  hist2[sout.str()]->SetLineColor(8);
	  
	  key_gluino1.push_back(sout.str());


	  sout.str("");
	  sout<<"glu2_"<<i;
	  hist2[sout.str()]=new TH2D(sout.str().c_str(),
				     sout.str().c_str(),
				     50, -3, 3,
				     50, 0, 3.142*2);
	  hist2[sout.str()]->SetOption("lego");
	  hist2[sout.str()]->SetLineColor(9);
	  key_gluino2.push_back(sout.str());
	}
    }
  //enable sum of squares weighing
  if(b_weight)
    {
      map<string, Plot>::iterator iter;
      for(iter = hist.begin(); iter!=hist.end(); ++iter)
	(iter->second)->Sumw2();

      map<string, Plot2>::iterator iter2;
      for(iter2 = hist2.begin(); iter2!=hist2.end(); ++iter2)
	(iter2->second)->Sumw2();
    }


}


